Positive-P phase space method simulation in superradiant emission from a cascade 

atomic ensemble 
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The superradiant emission properties from an atomic ensemble with cascade level configuration 
is numerically simulated. The correlated spontaneous emissions (signal then idler fields) are purely 
stochastic processes which are initiated by quantum fluctuations. We utilize the positive-P phase 
space method to investigate the dynamics of the atoms and counter-propagating emissions. The 
light field intensities are calculated, and the signal-idler correlation function is studied for different 
optical depths of the atomic ensemble. Shorter correlation time scale for a denser atomic ensemble 
implies a broader spectral window needed to store or retrieve the idler pulse. 

PACS numbers: 42.50.Lc, 42.50.Gy, 02.50.Ey, 02.60.Lj 



I. INTRODUCTION 

A quantum communication network based on the dis- 
tribution and sharing of entangled states is potentially 
secure to eavesdropping and is therefore of great prac- 
tical interest [H-Q- A protocol for the realization of 
such a long distance system, known as the quantum re- 
peater, was proposed by Briegel et al. 0,@. A quantum 
repeater based on the use of atomic ensembles as mem- 
ory elements, distributed over the network, was subse- 
quently suggested by Duan, Lukin, Cirac and Zoller Q. 
The storage of information in the atomic ensembles in- 
volves the Raman scattering of an incident light beam 
from ground state atoms with the emission of a signal 
photon. The photon is correlated with the creation of a 
phased, ground-state, coherent excitation of the atomic 
ensemble. The information may be retrieved by a re- 
verse Raman scattering process, sending the excitation 
back to the initial atomic ground state and generating an 
idlerphoton directionally correlated with the signal pho- 
ton [7Hl5j]. In the alkali gases, the signal and the idler 
field wavelengths are in the near-infrared spectral region. 
This presents a wavelength mismatch with telecommu- 
nication wavelength optical fiber, which has a transmis- 
sion window at longer wavelengths (1.1-1.6 urn). It is 
this mismatch that motivates the search for alternative 
processes that can generate telecom wavelength photons 
correlated with atomic spin waves (l6j . 

This motivates the research presented in this article 
where we study multi-level atomic schemes in which the 
transition between the excited states is resonant with a 
telecom wavelength light field [l6j]. The basic problem 
is to harness the absorption and the emission of telecom 
photons while preserving quantum correlations between 
the atoms, which store information and the photons that 
carry along the optical fiber channel of the network. 

It is not common to have a telecom ground state transi- 
tion in atomic gases except for rare earth elements [U Gil 
or in an erbium-doped crystal (19|. However, a telecom 
wavelength (signal) can be generated from transitions be- 
tween excited levels in the alkali metals [la, SOI ■ 



The ladder configuration of atomic levels provides 
a source for telecom photons (signal) from the upper 
atomic transition. For rubidium and cesium atoms, the 
signal field has the range around 1.3-1.5 /zm that can be 
coupled to an optical fiber and transmitted to a remote 
location. Cascade emission may result in pairs of pho- 
tons, the signal entangled with the subsequently emitted 
infrared photon (idler) from the lower atomic transition. 
Entangled signal and idler photons were generated from a 
phase-matched four-wave mixing configuration in a cold, 
optically thick 85 Rb ensemble [l(|. This correlated two- 
photon source is potentially useful as the signal field has 
telecom wavelength. 

The temporal emission characteristics of the idler field, 
generated on the lower arm of the cascade transition, 
were observed in measurements of the joint signal-idler 
correlation function. The idler decay time was shorter 
than the natural atomic decay time and dependent on 
optical thickness in a way reminiscent of superradiance 

011. 

The spontaneous emission from an optically dense 
atomic ensemble is a many-body problem due to the ra- 
diative coupling between atoms. This coupling is re- 
sponsible for the phenomenon of superradiance firstly 
discussed by Dicke HH in 1954. 

Since then, this collective emission has been exten- 
sively studied in two atom systems indicating a dipole- 
dipole interaction [U , in the totally inverted N atom 
systems (2ol \2l\ . and in the extended atomic ensemble 
[23j | . The emission intensity has been investigated using 
the master equation app roach [28U3Cj ] and with Maxwell- 
Bloch equations [U [32| . A useful summary and review 
of superradiance can be found in the reference (33l . [34| . 
Recent approaches to superradiance include the quantum 
trajectory method (35l . l36j and the quantum correction 
method [37| . 

In the limit of single atomic excitation, superradiant 
emission characteristics have been discussed in the refer- 
ence HH and [3{| . For a singly excited system, the basis 
set reduces to N rather than 2 N states. Radiative phe- 
nomena have been investigated using dynamical methods 
[40M42I ] and by the numerical solution of an eigenvalue 
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problem (43l - l46j . A collective frequency shift [13, EH 
can be significant at a high atomic density [4!| and has 
been observed recently in an experiment where atoms are 
resonant with a planar cavity |50|. 

To account for multiple atomic excitations in the 
signal-idler emission from a cascade atomic ensemble, the 
Schrodinger's equation approach becomes cumbersome. 
An alternative theory of c-number Langevin equations is 
suitable for solution by stochastic simulations. 

Langevin equations were initially derived to describe 
Brownian motion (5lj . A fluctuating force is used to 
represent the random impacts of the environment on the 
Brownian particle. A given realization of the Langevin 
equation involves a trajectory perturbed by the random 
force. Ensemble averaging such trajectories provides a 
natural and direct way to investigate the dynamics of the 
stochastic variables. 

An essential element in the stochastic simulations is 
a proper characterization of the Langevin noises. These 
represent the quantum fluctuations responsible for the 
initiation of the spontaneous emission from the inverted 
[32l I52II54I , or pumped atomic system [H, [56| as in our 
case. 

The positive-P phase space method [57T - I6H ] is employed 
to derive the Fokker-Planck equations that lead directly 
to the c-number Langevin equations. The classical noise 
correlation functions, equivalently diffusion coefficients, 
are alternatively confirmed by use of the Einstein rela- 
tions |6#j66j . The c-number Langevin equations corre- 
spond to Ito-type stochastic differential equations that 
may be simulated numerically. The noise correlations 
can be represented either by using a square (6?| or a 
non-square "square root" diffusion matrix |6lj |. The ap- 
proach enables us to calculate normally-ordered quan- 
tities, signal-idler field intensities, and the second-order 
correlation function. The numerical approach involves a 
semi-implicit difference algorithm and shooting method 
68] to integrate the stochastic " Maxwell- B loch" equa- 
tions. 

Recently a new positive-P phase space method involv- 
ing a stochastic gauge function [g9| has been developed. 
This approach has an improved treatment of sampling 
errors and boundary errors in the treatment of quantum 
anharmonic oscillators [7(J [zU ■ It has also been applied 
to a many-body system of bosons [72| and fermions [73[ . 
In this paper, we follow the traditional positive-P repre- 
sentation method [74j . 

The remainder of this paper is organized as follows. 
In section II, we show the formalism of positive P- 
representation, and demonstrate the stochastic differen- 
tial equations of cascade emission (signal and idler) from 
an atomic ensemble. In section III we solve numerically 
for the dynamics of the atoms and counter-propagating 
signal and idler fields in a positive P-representation. We 
present results of signal and idler field intensities, and 
the signal-idler second order correlation function for dif- 
ferent optical depths of the atomic ensemble. Section 
IV presents our discussions and conclusions. In the 



appendix, we show the details in the derivations of c- 
number Langevin equations that are the foundation for 
numerical approaches of the cascade emission. In Ap- 
pendix A, we formulate the Hamiltonian, and derive 
the Fokker-Planck equations by characteristic functions 
75] in positive P-representation. Then corresponding c- 
number Langevin equations are derived, and the noise 
correlations are found from the diffusion coefficients in 
Fokker-Planck equations as shown in Appendix B. 



II. THEORY OF CASCADE EMISSION 

The phase space methods [HH that mainly include P- 
, Q-, and Wigner (W) representations are techniques of 
using classical analogues to study quantum systems, es- 
pecially harmonic oscillators. The eigenstate of har- 
monic oscillator is a coherent state that provides the basis 
expansion to construct various representations. P and 
Q-representation are associated respectively with evalu- 
ations of normal and anti-normal order correlations of 
creation and destruction operators. W-representation 
is invented for the purpose of describing symmetrically 
ordered creation and destruction operators. Since P- 
representation describes normally ordered quantities that 
are relevant in experiments, we are interested in inves- 
tigating one class of generalized P-representations, the 
positive P-representation that has semi-definite property 
in the diffusion process, which is important in describing 
quantum noise systems. 

Positive-P representation [74| [76[ is an extension to 
Glauber-Sudarshan P-representation that uses coherent 
state (| a)) as a basis expansion of density operator 
p. In terms of diagonal coherent states with a quasi- 
probability distribution, P(a,a*), a density operator in 
P-representation is 

p = f \a)(a\P(a,a*)d 2 a, (1) 
Jd 

where D represents the integration domain. The nor- 
malization condition of p, which is Tr{p}= I, indicates 
the normalization for P as well, J D P(a, a*)d 2 a = 1. 

Positive P-representation uses a non-diagonal coher- 
ent state expansion and the density operator can be ex- 
pressed as 

P = f A(a,j3)P(a,l3)d(i(a,j3), (2) 

where 



d[i(a, p) = d 2 ad 2 (3 and A(a, 0) = , (3) 

(j3*\a) 

and (/?* \a) in non-diagonal projection operators, A(a, (3), 
makes sure of the normalization condition in distribution 
function, P(a, /3). 
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FIG. 1: Four-level atomic ensemble interacting with two driv- 
ing lasers (solid) with Rabi frequencies Q, a and Qb- Signal and 
idler fields are labelled by a 3 and Oj, respectively and Ai and 
A2 are one and two-photon laser detunings. 



Any normally ordered observable can be deduced from 
the distribution function P(a, f3) that 

((a t ) m a™) = / (3 m a n P(a,(3)dfi(a,(3). (4) 
Jd 

A characteristic function X P (A a , A^) (Fourier- 
transformed distribution function in Glaubcr-Sudarshan 
P-representation but now is extended into a larger 
dimension) can help formulate distribution function, 
which is 

Xp(X a ,Xp)= [ e^ a+iX ^P(a,(3)d^(a,(3). (5) 

JD 

It is calculated from a normally ordered exponential op- 
erator E(X), 

X P (A Q , Xp) = Tr{pE(\)}, E(X) = e iX ^e iX " a . (6) 

Then a Fokker-Planck equation can be derived from 
the time derivative of characteristic function, 

% = | T r{pi? ( A)}^Tr{gi?(A)} (7) 



by Liouville equations, 

In laser theory (75[, a P-representation method is ex- 
tended to describe atomic and atom-field interaction sys- 
tems. When a large number of atoms is considered, 
which is indeed the case of the actual laser, a macroscopic 
variable can be defined. Then a generalized Fokker- 
Planck equation can be derived from characteristic func- 
tions by neglecting higher order terms that are propor- 
tional to the inverse of number of atoms. It is similar 
to our case when we solve light-matter interactions in 
an atomic ensemble that the large number cuts off the 
higher order terms in characteristic functions. 

We consider N cold atoms that are initially prepared in 
the ground state interacting with four independent elec- 
tromagnetic fields. As shown in Fig[TJ two driving lasers 
(of Rabi frequencies il a and f2f,) excite a ladder configura- 
tion |0) — s- |1) — > |2). Two quantum fields, signal a s and 
idler Sj, are generated spontaneously. We note that the 
spontaneous emission from the cascade driving scheme 
is a stochastic process due to the quantum fluctuations, 
unlike the diamond configuration where quantum noise 
can be neglected (77117^1 . 

The complete derivation of the c-number Langevin 
equations for cascade emission from the four-level atomic 
ensemble is described in Appendix A and B. After set- 
ting up the Hamiltonian, we follow the standard pro- 
cedure to construct the characteristic functions |75ll in 
Appendix A using the positive-P representation |58j . In 
Appendix B.l, the Fokker-Planck equation is found by di- 
rectly Fourier transforming the characteristic functions, 
and making a 1 /N z expansion. 

Finally the Ito stochastic differential equations are 
written down from inspection of the first-order derivative 
(drift term) and second-order derivative (diffusion term) 
in the Fokker-Planck equation. The equations are then 
written in dimensionless form by introducing the Arecchi- 
Courtens cooperation units [79j in Appendix B.2. From 
Eq. (|B10I) and the field equations that follow, these c- 
number Langevin equations in a co-moving frame are, 
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where (I) stands for Ito type SDE. 7^ is the stochastic 
variable that corresponds to the atomic populations of 
state |i) when i = j and to atomic coherence when i ^ j, 
and Tij are c-number Langevin noises. The remaining 
equations of motion, which close the set, can be found 
by replacing the above classical variables, 7r*,, — > tt 
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n jk , (E+.)* -> e- v (E.y 



jk 



t 

jk' 

and 



Fj k . Note that the atomic populations satisfy 
njj. The superscripts, dagger (f) for atomic vari- 



ables and (— ) for field variables, denote the independent 
variables, which is a feature of the positive-P represen- 
tation: there are double dimension spaces for each vari- 
able. These variables are complex conjugate to each 
other when ensemble averages are taken, for example 

M = (4)* and ) = <£-.)• . The doubled 

spaces allow the variables to explore trajectories outside 
the classical phase space. 

Before going further to discuss the numerical solution 
of the SDE, we point out that the diffusion matrix el- 
ements have been computed using Fokkcr-Planck equa- 
tions and by the Einstein relations discussed in Appendix 
B.2. This provides the important check on the lengthy 
derivations of the diffusion matrix elements we need for 
the simulations. 

The next step is to find expressions for the Langevin 
noises in terms of a non-square matrix B [6ll [76j . The 
matrix B is used to construct the symmetric diffusion 
matrix D(a) = B(a)B T (a) for a Ito SDE, 



dx\ = Mtrft)dt+y^B ij (t,W l )dWl{t) (I) (10) 



where &dt = dW^(t) (Wiener process) and (^i(t)^j(f)} = 
5ijS(t — t'). Note that B — > BS, where S is an orthogo- 
nal matrix (SS T — I), leaves D unchanged, so B is not 
unique. We could also construct a square matrix rep- 
resentation B [m [H, [6?} ■ This involves a procedure of 
matrix decomposition into a product of lower and upper 
triangular matrix factors. A Cholesky decomposition 
can be used to determine the B matrix elements succes- 
sively row by row. The downside of this procedure is that 



the B matrix elements must be differentiated in convert- 
ing the Ito SDE to its equivalent Stratonovich form for 
numerical solution. 

The Stratonovich SDE is necessary for the stability 
and the convergence of semi-implicit methods. Be- 
cause of the analytic difficulties in transforming to the 
Stratonovich form, we use instead the non-square form 
of 

In this case a typical B matrix element is a sum of 
terms, each one of which is a product of the square root 
of a diffusion matrix element with a unit strength real 
(if the diffusion matrix element is diagonal) or complex 
(if the diffusion matrix element is off-diagonal) Gaussian 
unit white noise. It is straightforward to check that a B 
matrix constructed in this way reproduces the required 
diffusion matrix D — BB T . 



As pointed out in the reference 



the transverse 



dipole-dipole interaction can be neglected and nonparax- 
ial spontaneous decay rate can be accounted for by a 
single atom decay rate if the atomic density is not too 
high. We are interested here in conditions where the en- 
semble length L is significant and propagation effects arc 
non-negligible, and the average distance between atoms 
d = y/V/N is larger than the transition wavelength A. 
The length scales satisfy A < d -C L, and we consider a 
pencil-like cylindrical atomic ensemble. The paraxial or 
one-dimensional assumption for field propagation is then 
valid, and the transverse dipole-dipole interaction is not 
important for the atomic density we focus here. 

The theory of cascade emission presented here provides 
the solid ground for simulations of fluctuations that ini- 
tiate the radiation process in the atomic ensemble. A 
proper way of treating fluctuations or noise correlations 
and formulating SDE requires an Ito form that is derived 
from the Fokker-Planck equation. An alternative but 
more straightforward approach by making quantum to 
classical correspondence in the quantum Langevin equa- 
tion does not guarantee an Ito type SDE. That is the 
reason we take the route of Fokker-Planck equation, and 
the coupled equations of Eq.® are the main results in 
this section. 
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III. RESULTS FOR SIGNAL, IDLER 
INTENSITIES, AND THE SECOND-ORDER 
CORRELATION FUNCTION 

There are several possible ways to integrate the differ- 
ential equation numerically. Three main categories of al- 
gorithm used arc forward (explicit), backward (implicit), 
and mid-point (semi- implicit) methods [68] . The forward 
difference method, which Euler or Runge-Kutta methods 
utilizes, is not guaranteed to converge in stochastic inte- 
grations 1801 . There it is shown that the semi- implicit 
method [81j is more robust in Stratonovich type SDE 
simulations [82[ . More extensive studies of the stability 
and convergence of SDE can be found in the reference 
(83j . The Stratonovich type SDE equivalent to the Ito 
type equation (fTUl) . is 

dx\ = [Ai{t, 4) - \ Y^l ^)^J B ^(*' ^)\ dt 

j k 

+ ^Bij(t^ t )dWl (Stratonovich), (11) 
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FIG. 2: (Color online) Time-varying pump fields and time 
evolution of atomic populations. (Left) The first pump field 
Q a (dotted-red) is a square pulse of duration 50 ns and Qt 
is continuous wave (solid- blue) . (Right) The time evolu- 
tion of the real part of populations for three atomic levels 
<Tn = (Q13) (dash dotted-red), 022 = (012) (dotted- blue), 
(733 = (an) (solid-green) at z = 0, L, and almost vanishing 
imaginary parts for all three of them, indicate convergence 
of the ensemble averages. Note that these atomic populations 
are uniform as a function of z. 



which has the same diffusion terms Bij, but with modi- 
fied drift terms. This "correction" term arises from the 
different definitions of stochastic integral in the Ito and 
Stratonovich calculus. 

At the end of Appendix C.3, we derive the "correc- 
tion" terms noted above. We then have 19 classical vari- 
ables including atomic populations, coherences, and two 
counter-propagating cascade fields. With 64 diffusion 
matrix elements and an associated 117 random numbers 
required to represent the instantaneous Langevin noises, 
we are ready to solve the equations numerically using the 
robust midpoint difference method. 

The problem we encounter here involves counter- 
propagating field equations in the space dimension and 
initial value type atomic equations in the time dimension. 
The counter-propagating field equations have a boundary 
condition specified at each end of the medium. This is a 
two-point boundary value problem, and a numerical ap- 
proach to its solution, the shooting method [68|, is used 
here. 

Any normally-ordered quantity (Q) can be derived by 

ensemble averages that (Q) = Y^LiQi/ 'R where Qi is 
the result for each realization. 

In this section, we present the second-order correlation 
function of signal-idler fields, and their intensity profiles. 
We define the intensities of signal and idler fields by 

l s {t) = (E-{t)E+{t)) , /,(*) = (E-(t)E+(t)) , (12) 

respectively, and the second-order signal-idler correlation 
function 

G s ,i(t,r) = (E-(t)Er(t + r)E+(t + T)E+(t)) (13) 

where r is the delay time of the idler field with respect a 
reference time t of the signal field. Since the correlation 



function is not stationary [84j, we choose t as the time 
when G s ^i is at its maximum. 

We consider a cigar shaped 85 Rb ensemble of radius 
0.25 mm and L — 3 mm. The operating conditions of the 
pump lasers are (0 OJ f^, Ai, A 2 ) = (0.4, 1, 1, 0)703 where 
fl a is the peak value of a 50 ns square pulse, and fib is 
the Rabi frequency of a continuous wave laser. Four- wave 
mixing condition (Afc = 0) is assumed. The four atomic 
levels are chosen as (|0), |1), I 2 ), |3» = (|5S 1/2 ,F=3), 
|5P 3/2 ,F=4), |4D 5/2 ,F=5), |5P 3/a ,F=4». The natural 
decay rate for atomic transition |1) — > |0) or |3) — > |0) 
is 701 = 703 = 1/26 ns and they have a wavelength 780 
nm. For atomic transition |2) —> |1) or \2) — > |3) is 
712 = 732 = O.I56703 [H[ with a telecom wavelength 
1.53/im. The scale factor of the coupling constants for 
signal and idler transitions is g s /gi — 0.775. 

We have investigated six different atomic densities 
from a dilute ensemble with an optical density (opd) of 
0.01 to a opd = 8.71. In FigH d and H we take the 
atomic density p = 10 10 cm~ 3 (opd = 2.18) for exam- 
ple, and the grid sizes for dimensionless time Ai = 4 and 
space Az = 0.0007 are chosen. The convergence of the 
grid spacings is fixed in practice by convergence to the 
signal intensity profile with an estimated relative error 
less than 0.5%. 

The temporal profiles of the exciting lasers are shown 
in the left panel of FigfSJ The atomic density is chosen 
as p = 10 10 cm -3 , and the cooperation time T c is 0.35 ns. 
The right panel shows time evolution of atomic popula- 
tions for levels |1), |2), and |3) at z = 0, L, that are spa- 
tially uniform. The populations are found by ensemble 
averaging the complex stochastic population variables. 
The imaginary parts of the ensemble averages tend to 
zero as the ensemble size is increased, and this is a useful 
indicator of convergence. In this example, the ensem- 
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FIG. 3: (Color online) Spatial-temporal intensity profiles of 
counter-propagating signal and idler fields. (a) At z = 0, 
real (dashed-blue) and imaginary (solid-red) parts of signal 
intensity, (b) At z — L, real (dash dotted-blue) and imag- 
inary (solid-red) parts of idler intensity. (c) and (d) are 
spatial-temporal profiles for signal and idler intensities respec- 
tively. Both intensities are normalized by the peak value of 
signal intensity that is 7.56 x 10 -12 i? 2 . Note that the idler 
fluctuations and its non-vanishing imaginary part indicate a 
relatively slower convergence compared with the signal inten- 
sity. The ensemble size was 8xl0 5 , and the atomic density 
p = 10 10 cm -3 . 



ble size was 8xl0 5 . The small rise after the pump pulse 
il a is turned off is due to the modulation caused by the 
pump pulse fib, which has a generalized Rabi frequency 
\/ A'2 + 40?. This influences also the intensity profiles 
and the correlation functions. 

In FigOJ we show counter-propagating signal (—z) and 
idler (+z) field intensities at the respective ends of the 
atomic ensemble and their spatial-temporal profiles re- 
spectively. The plots show the real and imaginary 
parts of the observables, and both are normalized to 
the peak value of signal intensity. Note that the char- 
acteristic field strength in terms of natural decay rate 
of the idler transition (703) and dipole moment (di) is 
(di /h)E c w 36.3703. The fluctuation in the real idler field 
intensity at z — L and non- vanishing imaginary part in- 
dicates a slower convergence compared to the signal field 
that has an almost vanishing imaginary part. The slow 
convergence is a practical limitation of the method. 

In Fig |3] (a), we show a contour plot of the second-order 
correlation function G Sl i(t s , ti) where f, > t s . In Figure^] 
(b), a section is shown through t s w 75 ns where G s ^ is at 
its maximum. The approximately exponential decay of 
G Sj i is clearly superradiant qualitatively consistent with 
the reference (l6| . The non- vanishing imaginary part of 
Ggj calculated by ensemble averaging is also shown in 
(b) and indicates a reasonable convergence after 8xl0 5 
realizations. 

In TableHJ we display numerical parameters of our sim- 
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FIG. 4: (Color online) Second-order correlation function 
G s .i(~ts,ti). The 2-D contour plot of the real part of G 3 ,% with 
a causal cut-off at t s = ti is shown in (a). The plot (b) 
gives a cross-section at t a = t m ~ 75 ns, which is normalized 
to the maximum of the real part (dashed-blue) of G 3 .i. The 
imaginary part (solid-red) of G s ,i is nearly vanishing, and the 
number of realizations is 8xl0 5 for p = 10 10 cm -3 . 



ulations for six different atomic densities. The number 
of dimensions in space and time is M t x M z with grid 
sizes (At, Az) in terms of cooperation time (T c ), length 
(L c ). The superradiant time scale (TV) is found by fit- 
ting G S} i to an exponential function (e _t / T/ ), with 95% 
confidence range. 

TABLE I: Numerical simulation parameters for different 
atomic densities p. Corresponding optical depth (opd), time 
and space grids (Mt x M z ) with grid sizes (At,Az) in terms 
of cooperation time (T c ) and length (L c ), and the fitted char- 
acteristic time Tf for G a ,i (see text). 



P(cm" 3 ) 


opd 


M t x M z 


At(T c ), 
Az{L c ) 


Tc(ns), 
L c (m) 


fitted Tf 
(ns) 


5x10' 


0.01 


111 x 42 


0.3, 5x10"° 


4.89, 1.47 


25.9 


5xl0 8 


0.11 


101 x 44 


0.9, 1.5xl0~ 4 


1.55, 0.46 


24.6 


5xl0 9 


1.09 
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3.1 



In Fig J51 the characteristic time scale is plotted as a 
function of atomic density and the factor Nfi, and shows 
faster decay for optically denser atomic ensembles. We 
also plot the timescale T\ = 7^ / '(N fi+1) (ns) where p i s 
the geometrical constant for a cylindrical ensemble |23| . 
The natural decay time 7^ = 26 ns corresponds to the 
D2 line of 85 Rb. The error bar indicates the deviation 
due to the fitting range from the peak of G St i to approx- 
imately 25% and 5% of the peak value. The results of 
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FIG. 5: (Color online) Characteristic timescales, Tf and Ti 
vs atomic density p and the superradiant enhancement factor 
N/j,. Tf (dotted-blue) is the fitted characteristic timescale 
for G s j(t s — t m ,ti — t m + t) where t m is chosen at its max- 
imum, as in Figure [4] The error bars indicate the fitting 
uncertainties. As a comparison, Ti—'y^ / (N p + 1) (dashed- 
black) is plotted where 7^ = 26 ns is the natural decay time 
of D2 line of 85 Rb atom, and p is the geometrical constant for 
a cylindrical atomic ensemble. The number of realizations 
is 4xl0 5 for p = 5 x 10 7 , 5 x 10 s , 5 x 10 9 cm" 3 , 8x10 s for 
p = 10 10 , 2 x 10 10 cm -3 , and 16xl0 5 for p = 4 x 10 10 cm -3 . 



simulations are in good qualitative agreement with the 
timescale of T\ that can be regarded as a superradiant 
time constant of lower transition in a two-photon cascade 
[&5l l84j . Tf approaches independent atom behavior at 
lower densities, which indicates no collective behavior as 
expected. We note here that our simulations involve mul- 
tiple excitations within the pumping condition similar to 
the experimental parameters [16| . The small deviation of 
Tf and T x might be due to the multiple emissions consid- 
ered in our simulations other than a two-photon source. 
On the other hand the close asymptotic dependence of 
atomic density or optical depth in Tf and T\ indicates a 
strong correlation between signal and idler fields due to 
the four- wave mixing condition as required and crucial in 
experiment (l6| . 

For larger opd atomic ensembles, larger statistical en- 
sembles are necessary for numerical simulations to con- 
verge. The integration of 8xl0 5 realizations used in 
the case of p — 10 10 cm" 3 consumes about 14 days with 
Matlab's parallel computing toolbox (function " parfor") 
with a Dell precision workstation T7400 (64-bit Quad- 
Core Intel Xeon processors). 



IV. DISCUSSION AND CONCLUSION 

The cascade atomic system studied here provides a 
source of telecommunication photons that are crucial for 
long distance quantum communication. We may take 
advantage of such low loss transmission bandwidth in 
the DLCZ protocol for a quantum repeater. The per- 



formance of the protocol relies on the efficiency of gen- 
erating the cascade emission pair, which is better for a 
larger optical depth of the prepared atomic ensemble. For 
other applications in quantum information science such 
as quantum swapping and quantum teleportation, the 
frequency space correlations also influence their success 
rates [86| . To utilize and implement the cascade emission 
in quantum communication, we characterize the emission 
properties, especially the signal-idler correlation function 
and its dependence on optical depths. Its superradiant 
timescale indicates a broader spectral distribution which 
saturates the storage efficiency of idler pulse in an auxil- 
iary atomic ensemble [16| by means of EIT (electromag- 
netic induced transparency). Therefore our calculation 
provides the minimal spectral window (1/Tf) of EIT to 
efficiently store and retrieve the idler pulse. 

In summary, we have derived c-number Langevin equa- 
tions in the positive-P representation for the cascade 
signal-idler emission process in an atomic ensemble. The 
equations are solved numerically by a stable and conver- 
gent semi-implicit difference method, while the counter- 
propagating spatial evolution is solved by implementing 
the shooting method. We investigate six different atomic 
densities readily obtainable in a magneto-optical trap ex- 
periment. Signal and idler field intensities and their 
correlation function are calculated by ensemble averages. 
Vanishing of the unphysical imaginary parts within some 
tolerance is used as a guide to convergence. We find 
an enhanced characteristic time scale for idler emission 
in the second-order correlation functions from a dense 
atomic ensemble, qualitatively consistent with the su- 
perradiance timescales used in a cylindrical dense atomic 
ensemble [H,[2|. 
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Appendix A: Hamiltonian and Characteristic 
functions in Positive P-representation method 

The Hamiltonian H is in Schrodinger picture, and we 
separate it into two parts where Hq is the free Hamilto- 
nian of the atomic ensemble and one dimensional counter- 
propagating signal and idler fields, and Hi is the inter- 
action Hamiltonian of atoms interacting with two clas- 
sical fields and two quantum fields (signal and idler) as 
shown in Fig[TJ Dipolc approximation of — d ■ E and 
rotating wave approximation (RWA) have been made to 
these interactions. Using the standard quantization of 
electromagnetic field [57} . we have 
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3 M 



M 



i=l l=-M l=-M 
M 



l=-M 



i.r 



M 



ik a zi—iuj a t 



l=-M 



+ n^aWe- 1 ^'- 1 ^ + h.c. 



M 

h ^2 [g s V2M+ la%a s 



-ik s zi 



l=-M 



9l V2M + ia l \a hl e tk ^ + h.c. 



where a 



E 



E 



JV. 



(Al) 



(A2) 



fa(t)dio£ (k a ) / (2h) , and f a is slow varying temporal pro- 
file without spatial dependence (ensemble scale much less 
than pulse length). g s = d,23£(k s )/H, £{k) — y / hu/2e V 



blad form where we have for the four-level atomic system, 
(-) = 

M N z 

— I 2ct oi P a oi ~ a oi a oi P ~ P a oi a oi \ 

l=-M M 

_i_ ^12 roa.M»' „ajMt _ ^aMI ^aM , 
2 

^ 32 ^ 32 ^ ^ 32 " 32 



■ To ~ AM * AMT "AMT "AM ~U,£t~AtJl 
- [2(7^ ^^12 - 012 012 P - P CT 12 012 J 

32 " 32 



+ ^f[26& l p&&« - a^a^P ~ *&']}. (A4) 
The characteristic functions can be calculated as 



X - Tr{£(AM 

dx 
dt 



(A5) 



Tr{^(A)|} 



dt> s P' 
(A6) 



dt 



2A/+1 ' 



and z r 

of propagation that is equally split into 2M + 1 ele- 
ments. Commutation relations of field operators are 

..c r ik-n.(zi-z,i) 



m = -M, ...,M, and L is the length ^~di^p ~ M-^M ^)sp^ 



Qt )A ( Qt >L + ( Q t ) A-L + ( )sd' 



Tr{S(A)-[^,p]}, 
jat, 

Tr{£(A)l[tf L ,p]}, 
Tr{£(A)^[i^_ L ,p]}, 



(A7) 



2M+1 ' 



[ai,a],] = diii j and the matrix u;«/ = 

accounts for field propagation by coupling the local mode 
operators where k n = 2im/L. Note that the Rabi fre- 
quency is half of the standard definition. 



where Hq = Ha + Hl, Ha is the atomic free evolution 
Hamiltonian, Hl is the Hamiltonian for laser fields, and 
Ha-l — Hj. The detail of derivations in various charac- 
teristic functions can be found in laser theory [75j or the- 
ory of light- atom interactions in atomic ensembles (78|. 



The normally ordered exponential operator is chosen 



as 



Appendix B: Stochastic Differential Equation 



A distribution function can be found by Fourier trans- 
forming the characteristic functions, 



E(\)=Y[E l (X), 



E l (X) 



: e lA 19 CT 01 e * A 18°"l2 e lX n a 02 e lX 16 a 13 e 



! t i\\ sit 



iA 15*03p iA 



13" 11 
I ~l 



«A 12 ct 22 pjAiiO-33 p iA ! 10 a- 32 p iAg(To3 p iAg5-5 3 . 



iA 6*12p iA 5 ,J 01p iA 4 a s ,lp i A3a Si! iA 2 Sl , iAjOi,; 



(A3) 



Aside from the atom- field interaction ff = ^ [-Hj /°] , 
when dissipation from vacuum is considered (single 
atomic decay), we can express them in terms of a Lind- 



/(«) 



f 



-ia-X 



then 



If 



df 
dt 

dx _ 



at — d{i\\) 
the boundary terms, we have — 



(2?r) n 

9x 



X (A)dA 1 ...dA n , (Bl) 



e -iS.X^ dh dK (B2) 



use integration by parts and neglect 



a T / where a 

minus sign is from iXp. Correspondingly, if = e lA(3 , 



we have ^ 
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1. Fokker-Planck equation 



Let 

dt J ^ 



[CaSw + C ( a ] - L 5w + C%_ L 5 W + £ sp 6 u >}f, (B3) 

and we may neglect higher order derivatives (third order 
and higher) in various £'s. The validity of truncation to 
second order is due to the expansion in the small param- 
eter 1/N Z . 

If the Fokker-Planck equation is 



2L- -—a f-—A -r 92 92 

dt ~ da a -' dp pi + 2^ dad/3 + dpda 



)D a pf 
(B4) 

where A and D are drift and diffusion terms then we have 
a corresponding classical Langevin equation 



a 9 (z,t) = l-a l ^- ikiZ ^\a n {z,t) = ±-a l n , 
a 12 {z,t) = -^-a[ 2 ,a 13 (z,t) = jj-a l 13 , 



au(z,t) 



— t — i(w23+A2)i ik a zi—ik b zi—ikiZi 



(B6) 



where e lAkz = e ik a zi-ik b z l -ik,z l +ik s z l ^ We note that 



i^ujiva\ = c-^—a\, - i^^iya^ — -c-^-a[, (B7) 



and a l = N z — a l 13 — a l 12 — ct\ l7 which will be used in 
later coupled equations. Also for the field variables, 



da 

~dl 



d£ 
~dt 



A 



(B5) 



with a correlation function (r a Tp) = 6(t — t')D a p. So 
now we can derive the equations of motion according to 
various Cs, but we postpone them and derivations of dif- 
fusion coefficients after the scaling is made for a dimen- 
sionless form in the next subsection. The demonstration 
of various £'s can be found in laser theory [75] or theory 
of light-atom interactions in atomic ensembles [78j . 



2. Slowly varying envelopes and scaled equations of 
motion 



Here we introduce the slowly varying envelopes and 
define our cross-grained collective atomic and field ob- 
servables, then finally transform the equations in a di- 
mensionless form for later numerical simulations. Define 
slow varying observables that 

a 5 (z,t) ee — a \e- lk « Zl + l ^\ a e (z,t) = ^-e^i+iutt^ 
I 



a>7(z, i) = -j-^ a 7 e i 

a 8 ( Z ,t) EE l. a l se -^at+i^t+ik a z l -ik i Zl^ 



E~{z,t) 
E+(z,t) 



di/H 



V2M + la\e 



di/h 



V2M + la[e 



(B8) 



where we use the idler dipole moment in signal field scal- 
ing for the purpose of scale-free atomic equation of mo- 
tions, so we need to keep in mind that in calculating 
signal intensity or correlation function, an extra factor of 
(di/ds) 2 needs to be taken care of. 

We choose the central frequency of signal and idler 
as lj s = W23 + A2,Wi = UJ3 where Ai = u) a — lo\ and 
A2 = w a + cjt — cj 2 - With a scaling of Arecchi-Courtens 
cooperation length [79}, we set up the units of time, 
length, and field strength in the following, 



L r — cT r 



1 / dfnuji 



1 I 

T c di/H 



(B9) 



Now the slowly varying and dimensionless equations of 
motion with Langevin noises in Ito's form are 



— a 5 = (iAi ^-)a 5 + i^ a (a - a 13 ) + in* b a 7 - ia 16 Ef + J" 5 , 

— a 6 = i(A 2 - Ai + i22i_Zi)5 6 - iQ* a a 7 + zft b (5 13 - a 12 ) + ia 8 E+ e - lAkz + J" 6 , 
—a 7 = (iA 2 - -|)a 7 - ^a«6 + «^&3 5 + ia 9 E+e~ lAkz - ia w Ef + J7, 
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— ai3 = -701«13 + 7l2«12 + ^a«ig - «^a a 5 ~ ^f>ai 8 + zf2£a 6 + -^13, 

d 

— a 12 = -72«12 + «^bai 8 - itt* b a 6 + ia 14 E+ e ~ lAkz - ia 10 E~e tAkz + T 12 , 

at 
d 

—an = -703^11 + 732«i2 ~ ia 14 E+e~ zAkz + ia 10 E- e tAkz + ia 15 E+ - ia 9 E- + F lu 

f) I 

— a 8 = -(iAi + 7oi + 7o3 )a 8 - ifi*5 9 - zft b 5 14 + ^5 6 £ s - e ^A ' £^ + i5 19 S+ + JT 8 , 

<9 ~ 703- r> ~ i — c- iAfcz , v~ ~ \jp+ i t- 

— ag = -— ag - m a a8 + za7-fr s e + i{a - u\\)EJ + Sg, 

-5 14 - -(iA 2 + 703 2 72 )5i4 - tfija 8 + i(5 12 - 5n)£;e iAb + z5 17 £+ + J" 14 , (BIO) 



where 72 = 712 + 732, and field propagation equations are 



( 9 9 )F- 



-iAfcz 



- J4, 



ia 9 + 



(Bll) 



where other Langevin noises can be found by using the 
correspondence, for example, J 7 ^ <H- J19. 

Before we proceed to formulate the diffusion coeffi- 
cients, we need to be careful about the scaling factor 
for the transformation to continuous variables when nu- 
merical simulation is applied. Take for example, 



where t^Jt is a unit transformation factor from the 

\9i\ 

signal field strength to the idler one. For a recognizable 
format of the above equations used in the main context, 
we change the labels in the below, 



"5 7T01, a 6 «-» 7T12, a 7 «-» 7T 2, «8 ^ ^13, "g <H> 7T03, 
aiO <r+ 7r 32 , 5n O 7r 33 , 5l2 «-» 7T 2 2, 5l3 TTlli 

ai4 -o- 7r| 2 , 5i5 aie <-> 7r| 3 , an «• 7rJ 2 , 

5i 8 o tt} 2 , aig o 7rJ l5 



(B12) 



where 7Ty is the stochastic variable that corresponds to 
the atomic populations of state \i) when i = j and to 
atomic coherence when i ^ j. Note that the associated 
c-number Langevin noises are changed accordingly. 
The Langevin noises are defined as 



F 7 (z,t) 



l pi 

5_ e -ik a zi+iu a t J7 Q f^ z ^ — ±J5_ e ik b Zi+iui b t 



N z 



pi pi 

1 7 „-ik a z l +ik b z l +iu} b t+iu} a t JT^^ ^ — 13 



z ±V z 

F 14 (z,t) = —r[ A e- t( - U23+A2 '> t e tk ° z '- lkbZ, - tk * Zl 7 



T 9 (z,t) 

T 4 {z,t) 



^e-* k > z ^\T 12 (z,t) = ^., 



di/h 
_9j_ 

d t /h 



V2M + le- iuJst T l 4 , 
V2M + le iUit T[ 



(B13) 



{T e (z,t)H{z' ,t')) 

\ „ik b zi+iui b t „—ik a z,i+iuj a t' /pi pi' 

- Jp e e \ L 6 L 5 



N? 



,ik b zi+ioj b t e —ik a zi+iu a t^Q e ik a z t -iuj a t 



ig t y/2M + le ikiZl a l 10 a[]S(t - t')5 w 
! [i(fl a T c )a 6 + ia w {Ef/E c )] - t')T c 



5{z - z')L 



where we have used lim 



(B14) 



M->co 



™±l8 w =5{z-z>),2M + 



1 = jj-, and N c = ^j^- is the cooperation number. Then 
we have the dimensionless form of diffusion coefficients. 



D, 



6.5 



D 6y5 = [iQ. a a e + ia w E+] . 



5{t - t')S{z - z') (B15) 
(B16) 



The dimensionless diffusion coefficients Da are 
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(i) £> 5 ,5 = -i2VL a a 5 \ D 5 . 6 = i{Vt a a 6 + a 10 Ef); D 5J = -ifl a a 7 ; D 5 , 8 = i(fl a a$ + (an - otia)E+); 
D 5 ,9 = -i(0 5 9 + a 5 E^); D 5 . M = -ia 16 Ef ; D 5A3 = ia 16 Ef; D hM = -ia ls E+; D 5<w = 712S12; 

(ii) -De,6 = -i2Q b a 6 ; D 6}8 = -itt b a s ; D 6 ,w = -iQ b aio; Ms, 13 = —iO^ai + 7oi5 6 ; 
E>6,i6 = -ia 7 E~ + 7oi5i ; D 6 , 1S = 7oiai 2 ; 

(iii) I>7,8 = -ia & Ef; D 7 . 9 = -ia 7 Ef; 

(iv) X>s,9 = -ia s E^; D SjW = iQ b (a 12 - an); D 8 ,n = ifi^au; L> 8 ,i2 = -ififeSu; 

£>8,i3 = -ift* a otQ + iuivEf + 7oi5 8 ; £> 8 ,i6 = iaLi$E+ - ia 9 E^ + 701S11 + 732OK12; As.is = ianEf + 7oi5i4; 

(v) -D 9 ,g = -i2a 9 Ef; D 9i i = iaio-Ej*"; A),i5 = 732«i2; 

(vi) Dio,io = -^aioS+e^^* 2 ; Z?io,ii = «(^fcSi 6 - a 7 E~r) + 7 03 5 10 ; £10,13 = -«Oh5i 6 ; 
-Dio,i4 = i^bSis - ift* b ae + 7o3"i2; £>io,i9 = ia 6 E~; 

(vii) D u ,ii = iauE+e~ lAkz - ia w E~e tAkz + i5 15 E+ - ia 9 E~ + 732X12 + 703S11; 
£11,12 = ta 10 E-e lAkz - ta 14 Ete- tAkz ~ l32 S 12 ; 

(viii) Di2,i2 = i^b5i 8 - iQ* b a e ~ ia w E~e lAkz + ia 14: E+ e~ lAkz + 7 2 5i 2 ; -D1243 = -? : ^fc5i8 + i^ b a e - 7i2<3i 2 ; 
(ix)L»i 3j i 3 = if2 a 5ig - itt* a a 5 + i^b^is ~ ^6^6 + 7oi«i3 + 712S12; 

un „iAkz. r, Iff* I ■— iAkz miw 

(Xj-L'3,8 = 1 |2«a6e ; ^3,9 = "j 72 la 7 e ■ (^ 17 j 

\9i\ 



Before going further to set up the stochastic differ- 
ential equation in the next subsection, we remark on 
the alternative method to derive the diffusion coefficients 
from the Heisenberg-Langevin approach with Einstein re- 
lations [6J-|66|, and it provides the important check for 
Fokker-Planck equations. We note here that a symmet- 
ric property of the diffusion coefficients is within Fokker- 
Planck equation, whereas the quantum diffusion coeffi- 
cients in quantum Langevin equation do not have sym- 
metric property simply because the quantum operators 
do not necessarily commute with each other. 

3. Ito and Stratonovich stochastic differential 
equations 

The c-number Langevin equations derived from 
Fokker-Planck equations have a direct cor resp ondence 
to Ito- type stochastic differential equations [U [58[ . In 
stochastic simulations, it is important to find the expres- 
sions of Langevin noises from diffusion coefficients. 

For any symmetric diffusion matrix D(a), it can always 
be factorized into 

D(a) = B{a)B T (a) (B18) 

where B — > BS (an orthogonal matrixS" that SS T = I) 
preserves the diffusion matrix so B is not unique. The 
matrix B is in terms of the Langevin noises where £idt — 
dW\ (Wiener process) and (&(*)£,•(*')) = S lJ S(t-t') and 



I 

the £i below is just a random number in Gaussian distri- 
bution with zero mean and unit variance. 

In numerical simulation, we use the semi-implicit al- 
gorithm that guarantees the stability and convergence in 
the integration of stochastic differential equations. So a 
transformation from Ito to Stratonovich-type stochastic 
differential equation is necessary, 

dx\ = Ai(t, x\)dt + B ij(t> xl)dWi (Ito) (B19) 
i 

dx\ = [Ai(t, 4) - \ B ik& ^t)-^-B lk (t, 5|)]dt 

j k 

+ ^2B l3 (t,xi)dWi (Stratonovich) (B20) 

j 

where a correction in drift term appears due to the trans- 
formation. 

In the end we have the full equations with 19 variables 
in the positive-P representation, 64 diffusion matrix ele- 
ments, and 117 noise terms (random number generators). 
Nonvanishing corrections in drift terms are only for 0.5, 
a 6 , a 9 , 5io, S11, 5i2, 5i3, and they are itl a /2, iQ b , iEf, 
iE+/2, (-3703 +732)/4, -72/4, (-5701 + 712)/ 4 respec- 
tively. 

The Langevin noises can be formulated as a non-square 
form [6ll [73, and in numerical simulations, we have a 
factor ^ N 1 AtAz for Langevin noises J- and N ^ tAz for 
correction terms. 
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